﻿#include "discrete_face_utils.hxx"
#include "idx_mesh_utils.hxx"
#include <algorithm>
#include <chrono>
#include <iostream>
#include <set>
#include <unordered_map>
#include <unordered_set>

bndy_type BoundaryAnalyzer::GetSingularType(const SPApar_pos& p, const double tol) {
    if (info_.is_singular_ule && SamePointJudgeBy3D(p, { info_.u_le, p.v }, tol))
        return bndy_low_u;
    else if (info_.is_singular_uri && SamePointJudgeBy3D(p, { info_.u_ri, p.v }, tol))
        return bndy_high_u;
    else if (info_.is_singular_vle && SamePointJudgeBy3D(p, { p.u, info_.v_le }, tol))
        return bndy_low_v;
    else if (info_.is_singular_vri && SamePointJudgeBy3D(p, { p.u, info_.v_ri }, tol))
        return bndy_high_v;
    return bndy_unknown;
}

void BoundaryAnalyzer::GetFitParam(SPApar_pos& pos, const double tol) {
    if (is_greater_than(pos.u, info_.u_ri) || is_less_than(pos.u, info_.u_le)) {
        auto x = info_.surf.param_period_u();
        pos.u = resetToInterval(pos.u, info_.u_le, info_.u_ri, info_.u_ri - info_.u_le);
    }
    if (is_greater_than(pos.v, info_.v_ri) || is_less_than(pos.v, info_.v_le)) {
        pos.v = resetToInterval(pos.v, info_.v_le, info_.v_ri, info_.v_ri - info_.v_le);
    }
    bndy_type point_type = GetSingularType(pos);
    switch (point_type) {
    case bndy_low_u:
        pos.u = info_.u_le;
        break;
    case bndy_high_u:
        pos.u = info_.u_ri;
        break;
    case bndy_low_v:
        pos.v = info_.v_le;
        break;
    case bndy_high_v:
        pos.v = info_.v_ri;
        break;
    default:
        break;
    }
    GetFitParam(pos, info_.surf.eval_position(pos), tol);
}
void BoundaryAnalyzer::GetFitParamBy2D(SPApar_pos& pos, const double tol) {
    if (is_greater_than(pos.u, info_.u_ri) || is_less_than(pos.u, info_.u_le)) {
        auto x = info_.surf.param_period_u();
        pos.u = resetToInterval(pos.u, info_.u_le, info_.u_ri, info_.u_ri - info_.u_le);
    }
    if (is_greater_than(pos.v, info_.v_ri) || is_less_than(pos.v, info_.v_le)) {
        pos.v = resetToInterval(pos.v, info_.v_le, info_.v_ri, info_.v_ri - info_.v_le);
    }
    bndy_type point_type = GetSingularType(pos);
    bool u_judge = true, v_judge = true;
    SPApar_pos param;
    switch (point_type) {
    case bndy_low_u:
    case bndy_high_u:
        v_judge = false;
        break;
    case bndy_low_v:
    case bndy_high_v:
        u_judge = false;
        break;
    default:
        break;
    }
    if (u_judge) {
        param = info_.u_ri - pos.u < info_.u_tolerance ? SPApar_pos(info_.u_ri, pos.v) : SPApar_pos(info_.u_le, pos.v);
        if (same_par_pos(pos, param, tol)) {
            pos = param;
        }
    }
    if (v_judge) {
        param = info_.v_ri - pos.v < info_.v_tolerance ? SPApar_pos(pos.u, info_.v_ri) : SPApar_pos(pos.u, info_.v_le);
        if (same_par_pos(pos, param, tol)) {
            pos = param;
        }
    }
}
void BoundaryAnalyzer::GetFitParam(SPApar_pos& pos, const SPAposition& pos_3d, const double tol) {
    bndy_type point_type = GetSingularType(pos);
    SPApar_pos param;
    SPAposition new_pos_3d;
    bool u_judge = true, v_judge = true;
    switch (point_type) {
    case bndy_low_u:
    case bndy_high_u:
        v_judge = false;
        break;
    case bndy_low_v:
    case bndy_high_v:
        u_judge = false;
        break;
    default:
        break;
    }
    if (u_judge) {
        param = info_.u_ri - pos.u < info_.u_tolerance ? SPApar_pos(info_.u_ri, pos.v) : SPApar_pos(info_.u_le, pos.v);
        new_pos_3d = info_.surf.eval_position(param);
        if (same_point(pos_3d, new_pos_3d, tol)) {
            pos = param;
        }
    }
    if (v_judge) {
        param = info_.v_ri - pos.v < info_.v_tolerance ? SPApar_pos(pos.u, info_.v_ri) : SPApar_pos(pos.u, info_.v_le);
        new_pos_3d = info_.surf.eval_position(param);
        if (same_point(pos_3d, new_pos_3d, tol)) {
            pos = param;
        }
    }
}

bool BoundaryAnalyzer::SamePointJudgeBy3D(const SPApar_pos& a, const SPApar_pos& b, const double tol) {
    auto x = info_.surf.eval_position(a);
    auto y = info_.surf.eval_position(b);
    return same_point(x, y, tol);
}

bool BoundaryAnalyzer::IsSingular(const SPApar_pos& p, const double tol) {
    return GetSingularType(p, tol) != bndy_unknown;
}

template <ParamDirect X> bool BoundaryAnalyzer::IsOnBoundary(const SPApar_pos& p, const double tol) {
    if constexpr (X == ParamDirect::U) {
        return (IsOnBoundary(p, bndy_low_u, tol) || IsOnBoundary(p, bndy_high_u, tol));
    }
    else {
        return (IsOnBoundary(p, bndy_low_v, tol) || IsOnBoundary(p, bndy_high_v, tol));
    }
}

bool BoundaryAnalyzer::IsOnBoundary(const SPApar_pos& p, const double tol) {
    return (IsOnBoundary(p, bndy_low_u, tol) || IsOnBoundary(p, bndy_high_u, tol) || IsOnBoundary(p, bndy_low_v, tol) || IsOnBoundary(p, bndy_high_v, tol));
}

bool BoundaryAnalyzer::IsOnBoundary(const SPApar_pos& p, bndy_type type, const double tol) {
    SPApar_pos param = p;
    SPAposition p_3d_pos = info_.surf.eval_position(param);
    switch (type) {
    case bndy_low_u:
        param.u = info_.u_le;
        break;
    case bndy_high_u:
        param.u = info_.u_ri;
        break;
    case bndy_low_v:
        param.v = info_.v_le;
        break;
    case bndy_high_v:
        param.v = info_.v_ri;
        break;
    default:
        break;
    }
    SPAposition p_3d_pos_bndy = info_.surf.eval_position(param);
    return same_point(p_3d_pos, p_3d_pos_bndy, tol);
}


bndy_type BoundaryAnalyzer::GetBoundaryType(const SPApar_pos& p, const double tol) {
    if (fabs(p.u - info_.u_le) < info_.u_tolerance && IsOnBoundary(p, bndy_low_u, tol)) return bndy_low_u;
    if (fabs(p.u - info_.u_ri) < info_.u_tolerance && IsOnBoundary(p, bndy_high_u, tol)) return bndy_high_u;
    if (fabs(p.v - info_.v_le) < info_.v_tolerance && IsOnBoundary(p, bndy_low_v, tol)) return bndy_low_v;
    if (fabs(p.v - info_.v_ri) < info_.v_tolerance && IsOnBoundary(p, bndy_high_v, tol)) return bndy_high_v;
    return bndy_unknown;
}

SPApar_pos BoundaryAnalyzer::GetClosestPoint(const SPApar_pos& p, const SPApar_pos& guess) {
    SPApar_pos res = guess;
    if (IsSingular(p) || IsSingular(guess)) return res;
    if (info_.is_periodic_u && std::fabs(p.u - guess.u) > info_.u_tolerance) {
        double diff = p.u - guess.u;
        if (diff > 0) {
            res.u += info_.u_ri - info_.u_le;
        }
        else if (diff < 0) {
            res.u -= info_.u_ri - info_.u_le;
        }
    }
    if (info_.is_periodic_v && std::fabs(p.v - guess.v) > info_.v_tolerance) {
        double diff = p.v - guess.v;
        if (diff > 0) {
            res.v += info_.v_ri - info_.v_le;
        }
        else if (diff < 0) {
            res.v -= info_.v_ri - info_.v_le;
        }
    }
    return res;
}

template <ParamDirect X> Direct BoundaryAnalyzer::GetDirectRelateBoundary(SPApar_dir dir_2d) {
    // 获取两个参数值对应的边界方向
    SPApar_dir boundary_dir;
    if constexpr (X == ParamDirect::U) {
        boundary_dir = SPApar_dir{ 0, 1 };
    }
    else {
        boundary_dir = SPApar_dir{ 1, 0 };
    }
    double val = dir_2d * boundary_dir;

    double cos_theta = (dir_2d % boundary_dir);
    if (is_zero_mch(val)) {
        if (is_equal(cos_theta, 1)) {  // theta = 0
            return Direct::right;
        }
        else {
            return Direct::left;
        }
    }
    else if (val > 0) {
        return Direct::right;
    }
    else {
        return Direct::left;
    }
}

BoundaryAnalyzer::BoundaryAnalyzer(FACE* _face) : info_(_face) {
}

template bool BoundaryAnalyzer::IsOnBoundary<ParamDirect::U>(const SPApar_pos&, const double);
template bool BoundaryAnalyzer::IsOnBoundary<ParamDirect::V>(const SPApar_pos&, const double);
template Direct BoundaryAnalyzer::GetDirectRelateBoundary<ParamDirect::U>(SPApar_dir dir_2d);
template Direct BoundaryAnalyzer::GetDirectRelateBoundary<ParamDirect::V>(SPApar_dir dir_2d);

bool IsEdgeDegenerate(EDGE* ed) {
    bool ret = false;
    if (ed->start() == ed->end())  // same start and end vertex
    {
        if (!ed->geometry() ||                                                                // Null curve
            (ed->geometry()->equation().type() == 3) ||                                       // degenerate_curve_type
            (!ed->geometry()->equation().periodic() && !ed->geometry()->equation().closed())  // open curve
            )
            ret = true;
    }
    return ret;
};


bool operator==(point a, point b) {
    return is_equal(a.x, b.x) && is_equal(a.y, b.y);
}

bool operator!=(point a, point b) {
    return !(a == b);
}

FaceterFaceInfo::FaceterFaceInfo(FACE* face, double u_ri, double u_le, double v_ri, double v_le, double u_tol, double v_tol, bool is_periodic_u, bool is_periodic_v)
    : u_ri(u_ri), u_le(u_le), v_ri(v_ri), v_le(v_le), u_tolerance(u_tol), v_tolerance(v_tol), is_periodic_u(is_periodic_u), is_periodic_v(is_periodic_v), surf(face->geometry()->equation_for_update()) {
}

FaceterFaceInfo::FaceterFaceInfo(FACE* face) : face_ptr_(face), surf(face->geometry()->equation_for_update()) {
    SPAinterval u_interval = recover_line_endpoint(face->geometry()->equation().param_range_u());
    SPAinterval v_interval = recover_line_endpoint(face->geometry()->equation().param_range_v());
    u_le = u_interval.start_pt();
    u_ri = u_interval.end_pt();
    v_le = v_interval.start_pt();
    v_ri = v_interval.end_pt();
    is_periodic_u = face->geometry()->equation().periodic_u();
    is_periodic_v = face->geometry()->equation().periodic_v();
    u_tolerance = face->geometry()->equation().param_range_u().length() * (1 / 3.0);
    v_tolerance = face->geometry()->equation().param_range_v().length() * (1 / 3.0);
    is_singular_ule = face->geometry()->equation().singular_u(u_le);
    is_singular_uri = face->geometry()->equation().singular_u(u_ri);
    is_singular_vle = face->geometry()->equation().singular_v(v_le);
    is_singular_vri = face->geometry()->equation().singular_v(v_ri);
    if (surf.closed_u()) {
        if (surf.singular_u(v_le)) {
            is_singular_vle = true;
        }
        else {  // B_7  B14的样例需要通过下列判断才能确定奇点情况
            curve* vle = surf.v_param_line(v_le);
            SPAinterval interval = vle->param_range();
            if (is_equal(vle->length(interval.start_pt(), interval.end_pt()), 0.0)) {
                is_singular_vle = true;
            }
        }
        if (surf.singular_u(v_ri)) {
            is_singular_vri = true;
        }
        else {
            curve* vri = surf.v_param_line(v_ri);
            SPAinterval interval = vri->param_range();
            if (is_equal(vri->length(interval.start_pt(), interval.end_pt()), 0.0)) {
                is_singular_vri = true;
            }
        }
    }
    if (surf.closed_v()) {
        if (surf.singular_u(u_le)) {
            is_singular_ule = true;
        }
        else {
            curve* ule = surf.v_param_line(u_le);
            SPAinterval interval = ule->param_range();
            if (is_equal(ule->length(interval.start_pt(), interval.end_pt()), 0.0)) {
                is_singular_ule = true;
            }
        }
        if (surf.singular_u(u_ri)) {
            is_singular_uri = true;
        }
        else {
            curve* uri = surf.v_param_line(u_ri);
            SPAinterval interval = uri->param_range();
            if (is_equal(uri->length(interval.start_pt(), interval.end_pt()), 0.0)) {
                is_singular_uri = true;
            }
        }
    }
    sense = face->sense();
};

point operator-(point a, point b) {
    return { a.x - b.x, a.y - b.y };
}

const double GetX(const point& p) {
    return p.x;
}

const double GetY(const point& p) {
    return p.y;
}
double resetToInterval(double x, double l, double r, double period) {
    double intervalLength = r - l;
    double adjustedX = fmod(x - l, period);
    if (adjustedX < 0) {
        adjustedX += period;
    }
    return l + fmod(adjustedX, intervalLength);
};
bool SamePointJudgeBy3D(FACE* f, const SPApar_pos& a, const SPApar_pos& b, const double tol) {
    auto x = f->geometry()->equation().eval_position(a);
    auto y = f->geometry()->equation().eval_position(b);
    return same_point(x, y, tol);
}

void PostProcessTriangle(FACE* f, double nt, double st, CDT::Triangulation<double>* pcdt) {
    for (int i = 0; i < pcdt->triangles.size(); i++) {
        const CDT::Triangle& t = pcdt->triangles[i];
        SPApar_pos pl[3];  // 三角形节点的参数值列表
        for (int j = 0; j < 3; j++) {
            pl[j].u = (pcdt->vertices)[t.vertices[j]].x;
            pl[j].v = (pcdt->vertices)[t.vertices[j]].y;
        }
        for (auto& neibor_id : t.neighbors) {
            if (neibor_id > pcdt->triangles.size()) continue;
            const CDT::Triangle& neibor = pcdt->triangles[neibor_id];
            SPApar_pos npl[3];
            for (int j = 0; j < 3; j++) {
                npl[j].u = (pcdt->vertices)[neibor.vertices[j]].x;
                npl[j].v = (pcdt->vertices)[neibor.vertices[j]].y;
            }
            // 获取邻居三角形上与t对着的的顶点
            const SPApar_pos& nopsv = npl[CDT::opposedVertexInd(neibor.neighbors, i)];

            // 获取t上与邻居三角形对着的顶点
            CDT::Index id = CDT::opposedVertexInd(t.neighbors, neibor_id);
            const SPApar_pos& opsv = pl[id];
            // 获取两个公共点
            const SPApar_pos& v1 = pl[CDT::cw(id)];
            const SPApar_pos& v2 = pl[CDT::ccw(id)];

            // 获取参数域两个对角线上的
            const SPAposition& pos_center = f->geometry()->equation().eval_position({ (v1.u + v2.u) / 2, (v1.v + v2.v) / 2 });
            const SPAposition& pos_center_main = f->geometry()->equation().eval_position({ (nopsv.u + opsv.u) / 2, (nopsv.v + opsv.v) / 2 });
            // 判断两个三角形参数域能否翻转
            if (pos_center != pos_center_main) {
                continue;
            }
            SPAposition center(pos_center.x(), pos_center.y(), pos_center.z());
            const SPAposition& pos_v1 = f->geometry()->equation().eval_position(v1);
            const SPAposition& pos_v2 = f->geometry()->equation().eval_position(v2);
            const SPAposition& pos_opsv = f->geometry()->equation().eval_position(opsv);
            const SPAposition& pos_nopsv = f->geometry()->equation().eval_position(nopsv);
            // 当前为副对角线状态
            const SPAposition& sub_diag_center = { (pos_v1.x() + pos_v2.x()) / 2, (pos_v1.y() + pos_v2.y()) / 2, (pos_v1.z() + pos_v2.z()) / 2 };

            const SPAposition& diag_center = { (pos_opsv.x() + pos_nopsv.x()) / 2, (pos_opsv.y() + pos_nopsv.y()) / 2, (pos_opsv.z() + pos_nopsv.z()) / 2 };
            // 换对角线
            if ((center - sub_diag_center).len_sq() > (center - diag_center).len_sq()) {
                pcdt->flipEdge(i, neibor_id);
                break;
            };
        }
    }
}
bool checkTriangle(FACE* f, double nt, double st, SPApar_pos* p) {
    std::vector<SPAposition> posl(5, SPAposition(0, 0, 0));
    // 三个顶点的坐标数组，3 * 3，i为顶点下标
    for (int i = 0; i < 3; ++i) {
        posl[i] = f->geometry()->equation().eval_position(p[i]);
        posl[3].set_x(posl[3].x() + posl[i].x() / 3.0);
        posl[3].set_y(posl[3].y() + posl[i].y() / 3.0);
        posl[3].set_z(posl[3].z() + posl[i].z() / 3.0);
    }
    // 重心的坐标 posl[3]在三角形上，posl[4]在曲面上。
    // posl[3] = (posl[0] + posl[1] + posl[2]) * (1.0 / 3);  // 在三角形上的几何坐标

    SPApar_pos temp_2d = p[1] + (p[0] - p[1]) / 2;
    SPAposition temp_3d = f->geometry()->equation().eval_position(temp_2d);
    SPAposition temp_comp = posl[1] + (posl[0] - posl[1]) / 2;
    SPApar_pos centerPar((p[0].u + p[1].u + p[2].u) / 3.0, (p[0].v + p[1].v + p[2].v) / 3.0);
    for (int i = 0; i < 3; i++) {
        for (int j = i + 1; j < 3; j++) {
            temp_2d = p[i] + (p[j] - p[i]) / 2;
            temp_3d = f->geometry()->equation().eval_position(temp_2d);
            temp_comp = posl[i] + (posl[j] - posl[i]) / 2;
            if (st * st < (temp_comp - temp_3d).len_sq()) return false;
        }
    }

    double deltauv = SPAresabs * SPAresabs;

    SPAvector v0 = posl[1] - posl[0];
    SPAvector v1 = posl[2] - posl[1];
    SPAvector v2 = posl[0] - posl[2];
    bool b0 = v0.len_sq() < deltauv;
    bool b1 = v1.len_sq() < deltauv;
    bool b2 = v2.len_sq() < deltauv;
    if (b0 || b1 || b2) return true;

    posl[4] = f->geometry()->equation().eval_position(centerPar);  // 参数面的重心的几何坐标

    // 判断距离容差 / 表面容差
    if (st * st < (posl[4] - posl[3]).len_sq()) {
        return false;
    }

    // 判断角度容差/法向容差
    SPAunit_vector nor = f->geometry()->equation().eval_normal(centerPar);

    v0 = normalise(v0);
    // v1 /= v1.len();
    v2 = normalise(v2);
    // if(f->sense()) nor = -nor;
    // norl[0] = temp1.Cross(temp2);  // 三角形法向

    // norl[1] = Vector3(nor.x(), nor.y(), nor.z());  // 曲面重心处法向
    // norl[1] = norl[1] * (1.0 / norl[1].length());
    SPAvector normal_approx = v0 * (-v2);
    if (normal_approx.len_sq() < deltauv) return true;
    if (f->sense()) nor = -nor;
    normal_approx = normalise(normal_approx);
    double dot = normal_approx % nor;
    if (dot > 1.0)
        dot = 1.0 - SPAresabs / 2;
    else if (dot < -1.0)
        dot = -1.0 + SPAresabs / 2;
    double theta = 180.0 * acos(dot) / M_PI;
    if (theta > nt && theta < 180 - nt) {
        return false;
    }
    return true;
}
bool fixTriangle(FACE* f, double nt, double st, CDT::Triangulation<double>* pcdt, std::vector<CDT::V2d<double>>* ppoints) {
    bool res = true;

    /**
     * @todo: 合适的递归方式
     * 如果使用flag判断的效率比较低，会重复判断
     *
     * 先尝试使用只加点的方式，让CDT再进行重新连接。
     */
    for (CDT::Triangle t : pcdt->triangles) {
        SPApar_pos pl[3];  // 三角形节点的参数值列表
        for (int ind = 0; ind < 3; ++ind) {
            pl[ind].u = (pcdt->vertices)[t.vertices[ind]].x;
            pl[ind].v = (pcdt->vertices)[t.vertices[ind]].y;
        }
        if (checkTriangle(f, nt, st, pl)) {
            continue;
        }
        else {
            res = false;
            // 添加点.
            // 判断三角形的一边是不是在面的EDGE上
            SPApar_pos mid(0.0, 0.0);
            for (int i = 0; i < 3; ++i) {
                mid.u += pl[i].u / 3.0;
                mid.v += pl[i].v / 3.0;
            }
            ppoints->push_back(CDT::V2d<double>(mid.u, mid.v));
        }
    }
    return res;
}

void NewIndexedMesh(FACE* f, const std::vector<UV_VERTEX>& vertice, const CDT::TriangleVec& triangles) {
    MESH_MANAGER* mesh_manager;
    api_get_mesh_manager(mesh_manager);
    mesh_manager->begin_mesh_output(f, nullptr, nullptr);
    std::vector<void*> return_node_ptr(vertice.size());
    if (mesh_manager->need_counts()) mesh_manager->announce_counts(triangles.size(), vertice.size(), triangles.size() * 3);

    // 遍历所有顶点，添加到网格中
    const surface& suf = f->geometry()->equation();
    for (int i = 0; i < vertice.size(); i++) {
        const UV_VERTEX& pp = vertice[i];
        // 添加顶点到网格
        return_node_ptr[i] = mesh_manager->announce_indexed_node(i, pp.uv, pp.position, pp.normal);
    }
    // 遍历所有三角形，添加到网格中
    for (int i = 0; i < triangles.size(); i++) {
        int p1 = triangles[i].vertices[0];
        int p2 = triangles[i].vertices[1];
        int p3 = triangles[i].vertices[2];
        // 创建一个有三个顶点的polygon
        mesh_manager->start_indexed_polygon(i, 3);
        // 添加这个polygon对应的点
        mesh_manager->announce_indexed_polynode(i, 0, return_node_ptr[p1]);
        mesh_manager->announce_indexed_polynode(i, 1, return_node_ptr[p2]);
        mesh_manager->announce_indexed_polynode(i, 2, return_node_ptr[p3]);
        mesh_manager->end_indexed_polygon(i);
    }
    mesh_manager->end_mesh_output(f, nullptr, nullptr);
    mesh_manager->save_mesh_output(f, nullptr, nullptr);
}

void WritePointsToCSV(const std::string& filename, const point_vector& points) {
    std::ofstream file(filename);

    if (!file.is_open()) {
        std::cerr << "Failed to open file: " << filename << std::endl;
        return;
    }
    file << "index,x,y\n";
    file << std::setiosflags(std::ios::fixed) << std::setprecision(8);
    for (size_t i = 0; i < points.size(); i++) {
        file << i << "," << points[i].x << "," << points[i].y << "\n";
    }
    file.close();
}

void WriteTrianglesToCSV(const std::string& filename, const CDT::TriangleVec& triangles) {
    std::ofstream file(filename);

    if (!file.is_open()) {
        std::cerr << "Failed to open file: " << filename << std::endl;
        return;
    }
    file << "index,v0,v1,v2,n0,n1,n2\n";
    for (size_t i = 0; i < triangles.size(); i++) {
        file << i << "," << triangles[i].vertices[0] << "," << triangles[i].vertices[1] << "," << triangles[i].vertices[2] << "," << triangles[i].neighbors[0] << "," << triangles[i].neighbors[1] << "," << triangles[i].neighbors[2] << "\n";
    }
    file.close();
}

void SaveFaceToSAT(const char* filename, FACE* face) {
    FileInfo fileinfo;
    fileinfo.set_units(1.0);
    fileinfo.set_product_id("Example Application");
    api_set_file_info((FileIdent | FileUnits), fileinfo);
    api_set_int_option("sequence_save_files", 1);

    ENTITY_LIST entity_list_save;
    entity_list_save.add(face);

    FILE* file_ptr = fopen(filename, "w");
    if (file_ptr == nullptr) {
        std::cerr << "Failed to open file for writing: " << filename << std::endl;
    }

    // 设置保存选项
    AcisOptions* ao = NULL;    // 使用默认选项
    logical text_mode = TRUE;  // 以文本模式保存

    outcome result = api_save_entity_list(file_ptr, text_mode, entity_list_save, ao);
    if (!result.ok()) {
        std::cerr << "Failed to save entity list to file." << std::endl;
        fclose(file_ptr);
        return;
    }
    fclose(file_ptr);
}

bool CheckFitAndJudgeDiagonal(FACE* face, const Box& box, const double angle_tolerance, const double distance_tolerance) {
    // 根据主对角线或者副对角线构造三角形用于判断平面拟合
    const SPAposition& pos_botLeft = face->geometry()->equation().eval_position(box.getBotLeft());
    const SPAposition& pos_topRight = face->geometry()->equation().eval_position(box.getTopRight());
    const SPAposition& pos_botRight = face->geometry()->equation().eval_position(box.getBotRight());
    const SPAposition& pos_topLeft = face->geometry()->equation().eval_position(box.getTopLeft());
    const SPAposition& pos_center = face->geometry()->equation().eval_position(box.getCenter());

    // 三维空间上主对角线上的中点
    const SPAposition& main = pos_botRight + (pos_topLeft - pos_botRight) / 2;
    // 三维空间上副对角线上的中点
    const SPAposition& second = pos_botLeft + (pos_topRight - pos_botLeft) / 2;
    // 如果主对角线中点到模型的距离小于副对角线中点到模型的距离，则采用主对角线划分方式判断拟合
    SPApar_pos tmp_p0[3] = { box.getTopLeft(), box.getBotLeft(), box.getBotRight() };   // tmp_p0为左下三角形
    SPApar_pos tmp_p1[3] = { box.getTopLeft(), box.getTopRight(), box.getBotRight() };  // tmp_p1为右上三角形
    SPApar_pos tmp_p2[3] = { box.getTopLeft(), box.getTopRight(), box.getBotLeft() };   // tmp_p2为左上三角形
    SPApar_pos tmp_p3[3] = { box.getBotLeft(), box.getBotRight(), box.getTopRight() };  // tmp_p3为右下三角形
    if (checkTriangle(face, angle_tolerance, distance_tolerance, tmp_p0) && checkTriangle(face, angle_tolerance, distance_tolerance, tmp_p1) && checkTriangle(face, angle_tolerance, distance_tolerance, tmp_p2) &&
        checkTriangle(face, angle_tolerance, distance_tolerance, tmp_p3)) {
        return true;
    }
    return false;
}

void RemoveDuilplicatePointsAndEdges(point_vector& points, edge_vector& edges) {
    int idp = 0;                                          // 用于给不重复的点分配新的索引
    std::unordered_map<CDT::V2d<double>, int> mp_points;  // 点到新索引的映射
    std::vector<int> points_mapping(points.size(), 0);    // 原始点索引到新点索引的映射
    CDT::EdgeUSet mp_edge;                                // 用于存储不重复的边
    point_vector final_points;                            // 不重复的点集
    edge_vector final_edges;                              // 不重复的边集

    for (int i = 0; i < points.size(); i++) {
        if (mp_points.find(points[i]) == mp_points.end()) {  // 如果点不在映射中，则添加
            mp_points[points[i]] = idp++;
            final_points.push_back(points[i]);  // 添加到不重复的点集中
        }
        points_mapping[i] = mp_points[points[i]];  // 更新点的映射
    }
    for (int i = 0; i < edges.size(); i++) {
        auto [x, y] = edges[i].verts();                            // 获取边的两个顶点
        CDT::Edge new_edge(points_mapping[x], points_mapping[y]);  // 根据点的新索引创建新边
        if (mp_edge.find(new_edge) == mp_edge.end()) {              // 如果边不在集合中，则添加
            mp_edge.insert(new_edge);
            final_edges.push_back(new_edge);  // 添加到不重复的边集中
        }
    }
    points = std::move(final_points);  // 更新点集
    edges = std::move(final_edges);    // 更新边集
    return;
}
template <typename Fn_1, typename Fn_2>
std::unordered_map<CDT::Edge, std::vector<double>> CalculatePolygonEdgeIntersections(const point_vector& face_points, const edge_vector& face_edges, const edge_vector& polygon_edges, const point_vector& polygon_points, Fn_1 get_sweep_cord, Fn_2 get_cord) {
    using Index = int;      // 定义索引类型为整数
    using Coord = double;   // 定义坐标类型为双精度浮点数
    enum class EventType {  // 定义事件类型枚举类
        Start,              // 扫描开始事件
        End,                // 扫描结束事件
        Sweep               // 扫描线事件
    };
    struct EventNode {   // 定义事件节点结构
        Coord cord;      // 坐标
        Index edge_idx;  // 边的索引
        EventType type;  // 事件类型
    };

    std::vector<EventNode> scan_events;                                 // 扫描线事件列表
    std::unordered_map<Coord, std::vector<Index>> sweep_line_to_edges;  // 扫描线到边的映射

    for (int i = 0; i < face_edges.size(); i++) {      // 遍历所有边
        auto [idx_1, idx_2] = face_edges[i].verts();  // 获取边的两个顶点索引
        sweep_line_to_edges[get_sweep_cord(face_points[idx_1])].push_back(i);
    }

    for (auto& x : sweep_line_to_edges) {                                                                        // 对每个扫描线上的边进行排序
        std::sort(x.second.begin(), x.second.end(), [&](Index one, Index two) {                                // 根据边的另一个顶点的坐标排序
            return get_cord(face_points[face_edges[one].v1()]) < get_cord(face_points[face_edges[two].v1()]);  // 按坐标排序
            });
    }

    for (int i = 0; i < polygon_edges.size(); i++) {                                          // 遍历多边形的所有边
        auto [idx_1, idx_2] = polygon_edges[i].verts();                                      // 获取边的两个顶点索引
        if (get_sweep_cord(polygon_points[idx_1]) > get_sweep_cord(polygon_points[idx_2])) {  // 如果第一个点的坐标大于第二个点的
            std::swap(idx_1, idx_2);                                                         // 交换两个点
        }
        scan_events.emplace_back(get_sweep_cord(polygon_points[idx_1]), i, EventType::Start);  // 添加扫描开始事件
        scan_events.emplace_back(get_sweep_cord(polygon_points[idx_2]), i, EventType::End);    // 添加扫描结束事件
    }

    for (auto& sweep_line_cord : sweep_line_to_edges) {                           // 对每个扫描线坐标添加扫描事件
        scan_events.emplace_back(sweep_line_cord.first, -1, EventType::Sweep);  // 添加扫描线事件
    }

    std::unordered_map<CDT::Edge, std::vector<double>> result;                           // 结果映射，边到交点坐标的列表
    std::set<Index> active_edges;                                                        // 活动边集合
    std::sort(scan_events.begin(), scan_events.end(), [&](EventNode& a, EventNode& b) {  // 对所有事件按坐标排序
        if (is_equal(a.cord, b.cord)) {                                                   // 如果坐标相等
            return static_cast<int>(a.type) < static_cast<int>(b.type);                  // 按事件类型排序
        }
        return a.cord < b.cord;  // 按坐标排序
        });

    for (int i = 0; i < scan_events.size(); i++) {                // 遍历所有事件
        if (scan_events[i].type == EventType::Sweep) {            // 如果是扫描线事件
            Coord sweep_line_cord = scan_events[i].cord;         // 获取扫描线坐标
            for (auto x : active_edges) {                          // 遍历所有活动边
                auto [idx_1, idx_2] = polygon_edges[x].verts();  // 获取边的两个顶点索引
                Coord sweep_cord_1 = get_sweep_cord(polygon_points[idx_1]);
                Coord sweep_cord_2 = get_sweep_cord(polygon_points[idx_2]);
                Coord cord_1 = get_cord(polygon_points[idx_1]);
                Coord cord_2 = get_cord(polygon_points[idx_2]);
                if (is_equal(sweep_cord_1, sweep_cord_2)) continue;                                                                            // 如果两个顶点的sweep坐标相等，跳过
                double intersect_cord = cord_1 + (sweep_line_cord - sweep_cord_1) * (cord_2 - cord_1) / (sweep_cord_2 - sweep_cord_1);        // 计算交点坐标
                auto sweep_edges = sweep_line_to_edges[sweep_line_cord];                                                                      // 获取当前扫描线上的边
                for (auto edge : sweep_edges) {                                                                                                 // 遍历这些边
                    double low_cord = std::min(get_cord(face_points[face_edges[edge].v1()]), get_cord(face_points[face_edges[edge].v2()]));   // 获取边的最小坐标
                    double high_cord = std::max(get_cord(face_points[face_edges[edge].v1()]), get_cord(face_points[face_edges[edge].v2()]));  // 获取边的最大坐标
                    // 此处可以采用单调剖分 + 链表 优化
                    if ((is_equal(low_cord, intersect_cord) || is_less_than(low_cord, intersect_cord)) && (is_equal(high_cord, intersect_cord) || is_greater_than(high_cord, intersect_cord))) {  // 如果新坐标在边的坐标范围内
                        result[face_edges[edge]].push_back(intersect_cord);                                                                                                                      // 将新坐标添加到结果中
                    }
                }
            }
        }
        if (scan_events[i].type == EventType::Start) {      // 如果是扫描开始事件
            active_edges.insert(scan_events[i].edge_idx);  // 将边添加到活动边集合
        }
        if (scan_events[i].type == EventType::End) {       // 如果是扫描结束事件
            active_edges.erase(scan_events[i].edge_idx);  // 从活动边集合中移除边
        }
    }
    for (auto& x : result) {  // 对结果进行排序
        std::sort(x.second.begin(), x.second.end());
    }
    return result;  // 返回结果
}
template std::unordered_map<CDT::Edge, std::vector<double>> CalculatePolygonEdgeIntersections(const point_vector&, const edge_vector&, const edge_vector&, const point_vector&, decltype(GetX), decltype(GetY));
template std::unordered_map<CDT::Edge, std::vector<double>> CalculatePolygonEdgeIntersections(const point_vector&, const edge_vector&, const edge_vector&, const point_vector&, decltype(GetY), decltype(GetX));

template <typename Fn_1, typename Fn_2> std::unordered_set<int> JudgePointsInPolygon2D(const point_vector& points, const edge_vector& polygon_edges, const point_vector& polygon_points, Fn_1 get_sweep_cord, Fn_2 get_cord) {
    using Index = int;
    using Coord = double;
    enum class EventType {
        Start,  // polygon_edge的起点
        End,    // polygon_edge的终点
        Sweep   // 扫描线
    };
    struct EventNode {
        Coord cord;
        Index edge_idx;
        EventType type;
    };
    std::unordered_set<Coord> sweep_lines;                          // 扫描线
    std::unordered_map<Coord, std::vector<Coord>> intersect_cords;  // 用于存扫描线上产生的交点
    std::vector<EventNode> scan_events;                             // 扫描器
    for (const auto& p : points) {                                    // 获取每个点所在的扫描线
        sweep_lines.insert(get_sweep_cord(p));
    }
    for (int i = 0; i < polygon_edges.size(); i++) {
        auto [idx_1, idx_2] = polygon_edges[i].verts();
        if (get_sweep_cord(polygon_points[idx_1]) > get_sweep_cord(polygon_points[idx_2])) {
            std::swap(idx_1, idx_2);
        }
        scan_events.emplace_back(get_sweep_cord(polygon_points[idx_1]), i, EventType::Start);
        scan_events.emplace_back(get_sweep_cord(polygon_points[idx_2]), i, EventType::End);
    }
    for (auto& sweep_line_cord : sweep_lines) {
        scan_events.emplace_back(sweep_line_cord, -1, EventType::Sweep);  // 将每条扫描线抽象成一个点
    }
    std::sort(scan_events.begin(), scan_events.end(), [&](EventNode& a, EventNode& b) {
        if (is_equal(a.cord, b.cord)) {
            return static_cast<int>(a.type) < static_cast<int>(b.type);  // 保证扫描线正确的位置
        }
        return a.cord < b.cord;
        });                                      // 将所有的点按y进行排序
    std::unordered_set<Index> active_edges;  // 维护一个当前扫描线扫到的polygon
    for (int i = 0; i < scan_events.size(); i++) {
        if (scan_events[i].type == EventType::Start) {  // 如果当前的点为polygon的边的起点，则将该边加入到当前的扫描线中
            active_edges.insert(scan_events[i].edge_idx);
        }
        if (scan_events[i].type == EventType::Sweep) {  // 当前的点是扫描线，则计算当前扫描线与所有相交的边的交点
            Coord sweep_line_cord = scan_events[i].cord;
            for (auto x : active_edges) {
                auto [idx_1, idx_2] = polygon_edges[x].verts();  // 获取边的两个顶点索引
                Coord sweep_cord_1 = get_sweep_cord(polygon_points[idx_1]);
                Coord sweep_cord_2 = get_sweep_cord(polygon_points[idx_2]);
                Coord cord_1 = get_cord(polygon_points[idx_1]);
                Coord cord_2 = get_cord(polygon_points[idx_2]);
                if (is_equal(sweep_cord_1, sweep_cord_2)) continue;                                                                      // 如果两个顶点的sweep坐标相等，跳过
                double intersect_cord = cord_1 + (sweep_line_cord - sweep_cord_1) * (cord_2 - cord_1) / (sweep_cord_2 - sweep_cord_1);  // 计算交点坐标
                intersect_cords[sweep_line_cord].push_back(intersect_cord);                                                             // 由于扫描线的x或y值固定，只需要存y或x即可
            }
        }
        if (scan_events[i].type == EventType::End) {  // 如果当前的点为polygon的边的终点，则将该边从当前的扫描线中移除
            active_edges.erase(scan_events[i].edge_idx);
        }
    }
    std::unordered_set<Index> result_idx;  // 用于存储最终的结果
    for (auto& x : intersect_cords) {        // 对所有扫描线得到的交点进行排序 （这里可以用链表维护简单多边形的连续性，就可以不用对点进行排序了）
        std::sort(x.second.begin(), x.second.end());
    }
    for (int i = 0; i < points.size(); i++) {
        int cnt = 0;
        bool is_in_polygon = false;
        auto& cords = intersect_cords[get_sweep_cord(points[i])];  // 获取当前的点所在的扫描线的交点集合
        for (auto cord : cords) {
            if (is_equal(cord, get_cord(points[i]))) {
                is_in_polygon = true;
                break;
            }
            if (cord > get_cord(points[i])) {
                break;
            }
            cnt++;
        }
        if (is_in_polygon || !(cnt % 2)) {  // 根据cnt判断是否在面内
            result_idx.insert(i);
        }
    }
    return result_idx;
}
template std::unordered_set<int> JudgePointsInPolygon2D(const point_vector&, const edge_vector&, const point_vector&, decltype(GetX), decltype(GetY));
template std::unordered_set<int> JudgePointsInPolygon2D(const point_vector&, const edge_vector&, const point_vector&, decltype(GetY), decltype(GetX));

/// DEBUG
void DeBugOutputUVGraph::init(FACE* f) {
    SPAinterval u_interval = f->geometry()->equation().param_range_u();
    SPAinterval v_interval = f->geometry()->equation().param_range_v();
    b_u = std::make_pair(u_interval.start_pt(), u_interval.start_pt() + u_interval.length());
    b_v = std::make_pair(v_interval.start_pt(), v_interval.start_pt() + v_interval.length());
    dis_points.clear();
    triangulation_edges.clear();
    constraint_edges.clear();
    boundary_indices.clear();
    sample_indices.clear();
}
void DeBugOutputUVGraph::AddTriangles(point_vector points, CDT::EdgeUSet edges_set) {
    for (auto edge : edges_set) {
        dis_points.emplace_back(points[edge.v1()].x, points[edge.v1()].y);
        dis_points.emplace_back(points[edge.v2()].x, points[edge.v2()].y);
        triangulation_edges.emplace_back(dis_points.size() - 1, dis_points.size() - 2);
    }
}
void DeBugOutputUVGraph::AddconstraintEdge(std::pair<double, double> pointa, std::pair<double, double> pointb) {
    dis_points.push_back(pointa);
    dis_points.push_back(pointb);
    constraint_edges.emplace_back(dis_points.size() - 1, dis_points.size() - 2);
}
void DeBugOutputUVGraph::AddBoundaryPoint(std::pair<double, double> point) {
    dis_points.push_back(point);
    boundary_indices.push_back(dis_points.size());
}
void DeBugOutputUVGraph::AddSamplePoint(std::pair<double, double> point) {
    dis_points.push_back(point);
    sample_indices.push_back(dis_points.size());
}
void DeBugOutputUVGraph::output() {
    std::ofstream outfile("data.txt");
    outfile << "Uinterval:\n";
    outfile << b_u.first << " " << b_u.second << "\n";
    outfile << "Vinterval:\n";
    outfile << b_v.first << " " << b_v.second << "\n";

    // 输出二维点集数组
    outfile << "Points:\n";
    for (const auto& p : dis_points) {
        outfile << p.first << " " << p.second << "\n";
    }

    // 输出边界点索引数组
    outfile << "Boundary Indices:\n";
    for (int idx : boundary_indices) {
        outfile << idx << "\n";
    }

    // 输出边界点索引数组
    outfile << "Sample Indices:\n";
    for (int idx : sample_indices) {
        outfile << idx << "\n";
    }

    // 输出三角化后形成的边集合
    outfile << "Triangulation Edges:\n";
    for (const auto& e : triangulation_edges) {
        outfile << e.first << " " << e.second << "\n";
    }

    // 输出限制边集合
    outfile << "Constraint Edges:\n";
    for (const auto& e : constraint_edges) {
        outfile << e.first << " " << e.second << "\n";
    }

    outfile.close();
}

logical get_facet_edge_points_and_params(EDGE* edge, SPAposition*& pos_array, double*& param_array, int& num_pts) {
    AF_POINT* start;
    AF_POINT* end;
    AF_POINT* curr;
    int knt = 0;
    pos_array = NULL;
    param_array = NULL;

    // Find the list of AF_POINTS on the given edge, if one exists.
    if (AF_POINT::find(edge, 0, start, end)) {
        // Determine the number of points in the list.
        for (curr = start; curr != end; curr = curr->next(0)) {
            knt++;
        }
        knt++;

        // Allocate arrays of the proper size.
        pos_array = ACIS_NEW SPAposition[knt];
        param_array = ACIS_NEW double[knt];

        // Populate the arrays.
        int index = 0;
        for (curr = start; index < knt; curr = curr->next(0)) {
            pos_array[index] = curr->get_position();
            param_array[index] = curr->get_parameter();
            index++;
        }
    }
    num_pts = knt;
    return TRUE;

}


outcome api_get_facet_edge_points_test(EDGE* edge, SPAposition*& polyline, int& num_pts, AcisOptions* ao) {
    AF_POINT* p1 = nullptr, * p2 = nullptr, * p = nullptr;
    logical res = AF_POINT::find((ENTITY*)edge, FORWARD, p1, p2);
    num_pts = 0;
    polyline = nullptr;
    if (res) {
        num_pts++;
        p = p1;
        while (p != p2) {
            num_pts++;
            p = p->next(FORWARD);
        }
        polyline = ACIS_NEW SPAposition[num_pts];
        p = p1;
        for (int i = 0; i < num_pts; i++, p = p->next(FORWARD)) {
            polyline[i] = p->get_position();
        }
    }
    /** @todo: 补充返回信息 */
    return outcome();
}

SPAinterval recover_line_endpoint(const SPAinterval& inter) {
    const double INF = 1e9;
    if (inter.type() == interval_finite) {
        return inter;
    }
    else if (inter.type() == interval_finite_below) {
        return { inter.start_pt(), INF };
    }
    else if (inter.type() == interval_finite_above) {
        return { -INF, inter.end_pt() };
    }
    else if (inter.type() == interval_infinite) {
        return { -INF, INF };
    }
    return 0;
}


void get_face_tolerance(double in_body_diagnoal, facet_options_visualization* fo, double& out_face_distance_tolerance, double& out_face_angle_tolerance) {
    face_quality_level this_face_quality = fo->get_face_quality();
    if (this_face_quality == coarse) {
        out_face_distance_tolerance = in_body_diagnoal * 0.004;
        out_face_angle_tolerance = 40.0;
    }
    else if (this_face_quality == medium_coarse) {  // default
        out_face_distance_tolerance = in_body_diagnoal * 0.002;
        out_face_angle_tolerance = 30.0;
    }
    else if (this_face_quality == medium_fine) {
        out_face_distance_tolerance = in_body_diagnoal * 0.001;
        out_face_angle_tolerance = 20.0;
    }
    else if (this_face_quality == fine) {
        out_face_distance_tolerance = in_body_diagnoal * 0.0005;
        out_face_angle_tolerance = 10.0;
    }
}

void get_edge_tolerance(double in_body_diagnoal, facet_options_visualization* fo, double& out_edge_distance_tolerance, double& out_edge_angle_tolerance) {
    double out_face_distance_tolerance = 0, out_face_angle_tolerance = 0;
    get_face_tolerance(in_body_diagnoal, fo, out_face_distance_tolerance, out_face_angle_tolerance);
    edge_quality_level this_edge_quality = fo->get_edge_quality();
    if (this_edge_quality == medium) {
        out_edge_distance_tolerance = out_face_distance_tolerance;
        out_edge_angle_tolerance = out_face_angle_tolerance;
    }
    else if (this_edge_quality == better) {  // default
        out_edge_distance_tolerance = out_face_distance_tolerance / 2.0;
        out_edge_angle_tolerance = out_face_angle_tolerance / 2.0;
    }
    else if (this_edge_quality == best) {
        out_edge_distance_tolerance = out_face_distance_tolerance / 4.0;
        out_edge_angle_tolerance = out_face_angle_tolerance / 4.0;
    }
}